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Low Reynolds Number k-z and Empirical Transition Models 
for Oscillatory Pipe Flow and Heat Transfer 


Christopher Bauer 

ABSTRACT 


Stirling engine heat exchangers are shell-and-tube type with oscillatory flow (zero- 
mean velocity) for the inner fluid. This heat transfer process involves laminar-transition- 
turbulent flow motions under oscillatory flow conditions. A low Reynolds number k-z 
model, (Lam-Bremhorst form), was utilized in the present study to simulate fluid flow 
and heat transfer in a circular tube. An empirical transition model was used to activate 
the low Reynolds number k- z model at the appropriate time within the cycle for a given 
axial location within the tube. The computational results were compared with 
experimental flow and heat transfer data for; (1) velocity profiles, (2) kinetic energy of 
turbulence, (4) skin friction factor, (4) temperature profiles, and (5) wall heat flux. The 
experimental data were obtained for flow in a tube ( 38 mm diameter and 60 diameter 

long), with the maximum Reynolds number based on velocity being R^ mMX - 11840 , a 
dimensionless frequency (Valensi number) of Va = 80.2 , at three axial locations 
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XJD = 16 , 30 and 44 . The agreement between the computations and the experiment 
is excellent in the laminar portion of the cycle and good in the turbulent portion. 
Moreover, the location of transition was predicted accurately. The Low Reynolds 
Number k-z model, together with an empirical transition model, is proposed herein to 
generate the wall heat flux values at different operating parameters than the experimental 
conditions. Those computational data can be used for testing the much simpler and less 
accurate one dimensional models utilized in 1-D Stirling Engine design codes. 
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CHAPTER I 


Statement of the Problem 

(1.1) Introduction 

Since James Watt developed Steam Engine at the end of eighteenth century, 
people have been working very hard to develop new forms of engines to produce 
mechanical power. Starting from steam engines and numerous forms of internal 
combustion engines, presently researchers are working on developing the Stirling engine. 

The Stirling engine is a device whose operation is based upon the thermodynamic 
cycle proposed by Robert Stirling (1815), whose theoretical efficiency is the same as that 
off Carnot cycle. Its recent designs have been efficient, reliable, operating with a low 
noise level with a very long life. Furthermore, it can work in any atmospheric conditions 
and is pollution free. An interesting feature of the Stirling machine is that it can be 
designed to work as an electrical or mechanical power generator, or as a cooling machine 
(acting as a refrigeration device or as a heat pump for space heating). An important 
feature is that it can be driven by any heat source available in a particular application, 
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such as solar energy or even heat produced by the human body. 

There is a wide range of possible applications for Stirling machines, from a power 
source for future space stations, to a reliable motor for agricultural machinery, 
Stirling engine powered artificial heart. As an example, the view of one of the NASA’ 
research engines design for space applications is presented in Figure 1.1. 
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figure 1.1. Quarter Sectional View of NASA’s Stirling Space 


Engine SPDE. 


Power Development 


The 


most important aspect for the future development of Stirling engines is 


a very 





careful design of all parts where heat transfer process is taking place. Such design can 
be aided greatly by using computer simulation. The new computer programs presently 
used to design and model Stirling engines need improvement for application to wide 
variety of Stirling engines. Therefore, there is a strong need for better codes and more 
accurate information about flow field and heat transfer phenomena in the laminar, 
transition and turbulent regions of the flow. Because of the and complexity of the 
problem, which requires that much data be stored and processed, a CRAY supercomputer 
was used for the calculations. 

A schematic view of the Stirling engine, with locations of heat transfer and fluid 
flow problem areas identified, is presented in Figure 1.2. 
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Figure 1.2. Stirling Engine Schematic with Locations of Heat Transfer and Fluid Flow 
Problems Areas. 




A number of experiments have been performed at the University of Minnesota to 
understand laminar, transition and turbulent flows. In parallel with those experiments, 
jcveral numerical simulations of oscillatory flow with heat transfer have also been 
performed, (Patankar (1992); Koehler (1990); Ibrahim et. al. (1992); Kannapareddy 
( 1993 )). It is obvious that experiments are very difficult and expensive. Therefore, 
using a simulation program, a wider range of operating conditions in various geometrical 
configurations can be studied. 

Adding to the basic knowledge and understanding of transition and turbulent 
oscillating flow, present study can contribute to the development of a simulation program 
and the improvement of Stirling engine designs. 

(1.2) Literature Review 

Many studies discussing oscillating flow problems exist in the literature. In order 
to understand the character of this type of flow, previous experiments and numerical 
simulations will be discussed further here. 

Many different types of unsteady flow are characterized in the literature. Among 
them, are the periodic unsteady flows, such as pulsating flow and oscillating flow will 
be described here. These two types of cyclic flows can first be characterized by the 
value of the mean flow velocity. In oscillating flow, which can be a specific case of pul- 
°ting flow, the mean velocity within one cycle is zero. This means that the net mass 
trans ^ er within a cycle across any tube cross-section is also zero. In pulsating flow, these 
conditions are not necessarily satisfied. Pulsating flow can be defined as a superposition 


of a steady mean and an oscillatory flow. 

Richardson and Tyler (1929) first studied steady and oscillating flows. By 
examination of the sound waves in resonators, they found the velocity maxima not far 
away from the wall. 

Ohmi, et. al. (1982) examined a wide range from fully laminar to fully turbulent 
oscillating flows. They found that the velocity profiles in the laminar part of the cycle 
agree with the theoretical oscillating flow solution. But in the turbulent part of the cycle, 

they found that the velocity profiles agree well with the 1/7 power law for steady 
turbulent flow. 

Iguchi et. al. (1982) studied liquid oscillation in a U-tube. They determined that 
the change from laminar to transition flow takes place where the amplitude of oscillation 
begins to deviate from that predicted by the analytical solution for laminar flow. On the 
other hand, the change from transition to turbulent flow appears when the measured 
amplitudes of oscillation match those computed with the 1/7 power law profile. 
However, it is difficult to directly compare their prediction with the transition regime 
location in a straight tube, because they used U-type bend tube, rather then straight tube. 

Hino, et. al. (1983) conducted an extensive study of oscillating flow in a 
•^^tangular duct. They focused on the turbulence structure of the flow, wall shear 
$tresses > Reynolds shear stresses, turbulent fluctuations, and coherent structure of 
turbulence. They found that in the accelerating phase turbulence is triggered near the 

*311 but suppressed, while in the deceleration phase, turbulence is generated vigorously 
*n the near wall 


region and spreads to the core flow. 


In his experimental work, Seume (1988) defines the criterion for transition as a 
„pid increase in the measured rms velocity fluctuations. He verifies that in oscillating 
flow, the critical Reynolds number depends on the Valensi number, and describes two 
mechanisms that trigger turbulence. First, transition can be cause by an incoming fluid 

canying a turbulent slug, and second, can be triggered by unstable boundary layer 
growth at higher Reynolds numbers. 

Koehler (1990) used the low Reynolds number k-e turbulence model to 

numerically predict the oscillating flow and associated heat transfer. He identified the 

Lam-Bremhorst form of the k-t model (Lam and Bremhorst (1981)), as suitable for 

oscillating flow calculations and compared mean velocity profiles and fluctuations with 

experimental results obtained from the oscillating flow test facility at the University of 

Minnesota. He showed that the model has the capability to qualitatively correctly predict 

the transition to turbulence in quasi-steady and accelerated pipe flow . He pointed out 

that the inflow boundary condition should be theoretically, or, if possible, experimentally 

investigated in order to enable the prediction of traveling a turbulent slug downstream 
of the flow. 

Eckmann and Grotberg (1991) studied the transition to turbulence in oscillating 
flow a pipe over a wide range of relative amplitude of fluid displacement (d,) and 

ley parameter The y USK * two measurement techniques: hot-film 

"omometry and laser Doppler velocimetry. n,ey observed that post transition 
•“Aulence was confined to a thin region near the wall for Reynolds numbers (based on 
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the Stokes-layer thickness Re & = U mnm 6/v ) greater than 500 and high frequencies. 

Ahn and Ibrahim (1992) used the high Reynolds number k-e turbulence model 
to examine oscillating flows under an operating conditions typical of Stirling power 
systems. Their results were compared with experimental data from University of 
Minnesota. In the laminar flow regime their predictions matched the data with relatively 
high accuracy; in the transition and turbulent regimes the computations matched the data 
with acceptable error. However, they concluded that further improvement in the 
turbulence modeling was necessary. 

Ibrahim and al. (1992) proposed an empirical model for transition to turbulence 
in oscillatory flows, in straight tubes. They used the momentum thickness Reynolds 

number (Re 4j ) at the point of transition to turbulence expressed in terms of turbulent 

intensity (77). From that model, the position of the turbulent slug from the tube entrance 
could be determined for different angular positions within the cycle. 

In experiments at the University of Minnesota documented by Seume et al. 
(1992), oscillating flow study at the SPDE Stirling engine operating conditions were 
conducted. They measured and documented the axial and radial components of ensemble 
averaged velocity, rms velocity fluctuation, and dominant Reynolds shear stress, at 
various radial position for four axial locations. From this data, the laminar, the 
fransition and the turbulent regions within the cycle could be identified and isolated. 
Their detailed measurements are useful in characterizing attributes of oscillating flow, 
ujcluding flow phenomena observed in the near wall region at flow reversal and during 
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ihe transition process. 


(1.3) Outline of Present Study 

First, the problem to be studied is mathematically described. All governing 
equations, the boundary conditions and important assumptions are introduced. Next, the 
numerical methods for solving the system of governing equations are presented. Then, 
the turbulence modeling is discussed. Following this, in order to validate code 
performance the steady state flow calculations are given. Finally, the results of the 
oscillating fluid flow and associated heat transfer are compared with experimental data 
and discussed. 
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CHAPTER II 

Mathematical Description of the Physical Phenomenon 

The description of the governing equations for all dependent variables, the 
fundamental assumptions and the boundary conditions for solving the oscillatory flow 
with heat nansfer problem is given in this chapter. Also, several nondimensional 
cables which not only simplify the problem but also provide the natural scale for the 
boundary condrdons, phystcal properties, and geometry are presented. 

(2.1) Governing Equation 

The governing differential equation expressing the conservation of mass, 
momentum and energy for a continuum are listed below. In k-t turbulence modeling, 
additional equations for turbulent kinetic energy (* -equation) and dissipation rate 
^ turbulence (« -equation) are used. All of these are given in the axisymmetric 
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coordinate system used for unsteady flow over an infinitesimal control volume, (see Peric 
and Scheuerer (1989) and Munson et. al. (1990)). 


Continuity Equation 

From the conservation of mass principle the following formulation for the 
continuity equations results. 


dU_ 

dx 


r dr 


0 


( 2 . 1 ) 


Momentum Equations 

The conservation of momentum yields two equations in the x and r directions 
respectively. 


x -Momentum 
dt dx dr 




t -Momentum 
pf-p 

dt dx dr 


dP* . d , dV. Id, dV . V 

dr dx^* dx ) *^1*^*1^ 


The modified pressure P* is equal to 




(2.3) 


(2.4) 


P is the hydrodynamic pressure. 


£ nerev Equation 




where S T is a source or sink term, representing e.g. heat generation by chemical 
reactions. In present study S T is assumed to be zero. 


Turbulent Kinetic Energy k -Equation 

dk j . dk . -dk 9 f. dk ] 1 9 f . dk l 

■ &l (,1 ^ ) &] + 7¥[ (,1 ^ ) ¥]^ p<? - pe (2 ' 6) 


where k = (l/2)(I/ /2 + V* 1 * W' 2 ) 

Dissipation Rate of Turbulence e -Equation 

( 2 . 7 ) 


where the generation term G is 


or dr r dr dx 


( 2 . 8 ) 


*nd using Prandtl Komogorov’s expression, the turbulent viscosity is modeled as: 

k 2 


P r s pc / — 


( 2 . 9 ) 


* -equation and e -equation with all constants and functions will be fully described 
in Chapter IV. 


V-i) Basic Assumptions 

' The geometry of Stirling engine heat exchangers which, 

..... ® 316 normally of the shell 

and-tube type, is complicated. However r • . me shell- 

“ ° f " e “ . a tTetri thc 

University of Mmnesota asweiiasin the computations, phase of fte ^ J 

" asymmetric coordinate system f.,*« 8 

y <X ’ r ’*>’ on = can assume that there are no c 
in the azimuthal direction anrf ,t changes 

’ a " d vel ° d * ta Ws direction is zero: 


w 


0 ’ d<p = 0 


Throughout this work the fh.iw • 

’ We fluid 15 ^sumed to be Newtonian , 

continuum, and. With that i„ mi „d and negiectina bu ’ ' nC ° mpre “ lb,e ' "" » 

“ination ufcs the following form: * ^ 


dU 


P~^-pO(VU) = - V /> + 


v (n(vt7)) + V(p(vi7)*) 


WhCre 11,6 P ress ure is defined as: 
u »ng Fourier’s law 


P = P + jF(Vt7) 


( 2 . 11 ) 


? = -A(vr) 

illr ^ " inC0 “- ~ W— .as. can be wri,te n as ’ 
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( 2 . 13 ) 


the viscous dissipation function <& for a Newtonian fluid in axsisymmetric 
coordinates can be written as: 

* . (2.14) 

dx dr r or ax 


The equations (2.1), (2.11), and (2.13) are sufficient to solve for the four 
independent variables: U, V, T and P. These form the complete set of equations 
necessary to describe the flow field and heat transfer in oscillating flow. The basic 
assumptions used to derive the theoretical equations are summarized as foUows: 

• Axisymmetric geometry. 

• Fluid is a continuum, Newtonian fluid with constant properties. 

• Fourier heat conduction law, no internal heat sources and no radiation heat 
transfer. 

• No body forces or gravitational effects. 


(2.3) Dimensional Analysis (Similarity Parameters) 

The following characteristic parameters are generally used to describe the fluid 
flow and heat transfer in similar systems. 


Maximum Reynolds Number (-Rg m . T ) 

One of the most important similarity parameters is a maximum Reynolds number, 


Re = *' v mu^ 

This definition is similar in structure to the well kn 

*- »• mean velocity u u ^ defi " , ‘ , ° n ° fRey " 0,ds Number 

for oscillating flow th 

(amplitude velocity) „ typi^ emp , oy£d ’ # "■«% ^ 

Vaiensi Numh < » r ^ 

sciliatory flows are unsteady bv drfi • • 

°e described: 

Fa = ^-P 2 

The VaJetisi Number can be physically • 

^ inter Preted as the ratin e 

**** *’/4, to the oscillation . 

osculation penod J/ Q 

S! «BUSi !mfc£t 

The Srouhal Number in- ; 

■Sm is used to connect the r 

^ the frequency of the oscillating flow 


caj 


*Hh £/• 0 ■ 

*“ ^ aod is defined as follows: 


Sir = 


^ earlier definitions of * 

mu 2nd Fl 2 J talfe« ih r 

taJces die form: Str = 4 Fa 

ttUx 

IadVC am P Ji ^ e of fluid displacement ■ 

ISa drived naramoh • 

^< | ^i ,• Parameter that ca 
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be defined as the maximum axial fluid displacement during one cycle divided by the duct 
length. In the mathematical form, 



Using ear lier definitions of P^ mMX and Vu , it takes the form. 

A - 1 
r " 2 L Va 

For A r > 1 all of the fluid initially contained in the pipe is pushed outside at 

some time during the cycle; 

For A r = 1 all of the fluid initially contained in the pipe is moved the length 

of the channel; 

For A r < 1 from all of the fluid initially contained in the pipe some does not 

leave for any time during the cycle. 

Prandtl Number ( Pr ) 

The well known definition of the Prandtl Number Pr, is: 


Pr 



where X is used for thermal conductivity. 


Geometric Similarity Parameters (XI D) and (D) 


The similarity parameters XI D and D are used to describe the geometry of 
tbc problem. The dimensionless length XI D , is a well known similarity parameter and 
*111 be often used in this work. 


(2.4) Boundary Conditions 

Since the governing momentum and energy equations are parabolic in time and 
elliptic in space, they can describe a whole class of fluid flow problems. By specifying 
boundary conditions for each of the dependent variables along the computational domain, 
ihe problem formulation may be completed. Four different types of boundary condition 
are used herein. They are shown in Figure 2.1. 


symmetry line 



Figure 2.1. Boundary Conditions. 
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Solid Walls 

Along the wall the fluid velocity is zero, assuming a no-slip condition and an 
impermeable wall. For the energy equation the wall temperature is specified. 

’lr 0 and r=r -'f = 0 e- 15 ) 

Axis of Symmetry 

At the center line of the tube, which is an axis of symmetry Neumann type 
boundary conditions are used here. 

Inlet Plane 

Dirichlet type boundary conditions are used on the inlet plane since all values of 
the dependent variables are specified. . 

U * * V^inB , V m = 0 , r* = (2.17) 

Outlet Planfr 

Neumann type boundary conditions are used on the outlet plane, which states in 

effect that the diffusive fluxes normal to the exit plane can be neglected. Since the values 

°f the dependent variables are not known a priori at these planes the gradients of the 

dependent variables in a direction normal to the outlet plane are assumed to be zero. 

** a situation is valid if the outlet plane is far from the entrance or any recirculating 
activities. 

dU _ dV _ dT Q 

dx etc etc 


( 2 . 18 ) 


It should be noted that by specifying the inlet velocity IA, as a time dependent 

value, (see equation 2.17) the acceleration and deceleration of the fluid is accounted. 
Simulating the oscillating flows after each half cycle, the outlet plane boundary 
conditions are used for the inlet plane boundary conditions during the neat half cycle. 
The boundary conditions for additional differential equations are presented in Chapter IV. 
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CHAPTER III 


Numerical Methods for the Solution of the Governing Equations 

The governing elliptic partial differential equations together with the boundary 
conditions can not be solved analytically. However, solutions for this coupled set of 
equations can be found using numerical methods. Any of the numerical methods of 
solution consist of two basic steps. First, the computational domain is divided into a 
number of subdomains. This transforms (discretizes) the partial differential equations into 
an easy-to- solve set of algebraic equations. Second, this set of algebraic equations has 
to be solved. 

Many methods to discretize partial differential equations have been proposed. 
Among them, the most widely used are the Finite Difference Method, Finite Element 
Method and Finite Volume Method. 

The code used in present study is based on the C.A.S.T. (Computer Aided 
Simulation of Turbulent Flows) code developed by Peric and Scheuerer (1989). The 
Finite Volume Method is used in this code with an ability to handle various flow types, 
isometries and boundary conditions. 
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(3.1) Discretization Method 

An important step in solving a set of partial differential equations is transforming 
fccm into one general transport equation. Next, one common algorithm is needed to solve 
ill dependent variables. Then, we can write for any scalar variable <j> the transport 
equation as: 


d(p4>) , d 

dt dx 


|(pt'*-r ( %.ii( P rF*-rA.s ( (3.i) 

dx 9 dx r dr * dr 9 


Choosing appropriate quantities for <f> , and from table I, all the transport 
partial differential equations can be retrieved from equation (3.1). 


Equation 


Continuity 


Momentum 


Momentum 


Energy 


Turbulent Kinetic 
Energy 


Turbulent 
Dissipation Rate 



Cj/iPG— c 2 f 2 p— 


TaMe /. Interpretations of <J> , T + and for the governing equations. 


4 i 

















ftoc computational domain is divided into a number of finite control volumes which 
compose the computational grid. To prevent the probable development of wavelike 
velocity or pressure fields, a "staggered grid" for velocity U and V has been used. For 
*11 other variables, a "normal grid" has been used (see figure 3.1). The values of the 
dependent variables are evaluated at the centers of their control volumes. 



Figure 3.1. The staggered grid for three distinctive spatial control volumes for: 

a) x-momentum equation; 

b) r-momentum equation; 

c) continuity, energy, k-t equations. 

Simplified diagrams of a computational domain and grid arrangement, are presented on 
figures 3.2 and 3.3. For a more detailed explanation refer to Peric and Scheuerer (1989). 
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(3.2) Solution Procedure 

The dependent variables are coupled together in the set of governing equations 
presented in Chapter H. The iterative SIMPLE algorithm (Semi Implicit Method for 
Pressure Link Equations) developed by Patankar and Spalding (1972) has been 
implemented in the C.A.S.T. code (Peric and Scheuerer (1989)), and has been used in 
this study. The SIMPLE algorithm consists of the following steps: 

^ tep q) The initialization of all the values for dependent variables are performed 

in order to evaluate the finite volume coefficients. For unsteady flows, 
values from the previous time step can be used or the given initial values 
for the first time step. 

(step 1) The linear equation set is created by assembling the finite volume 

coefficients of the x -momentum equation. The resulting set is solved 

giving the velocity field U m . 

(step 2) The same procedure is performed for the r -momentum equation to 

give V . 

(step 3^ As the initially guessed pressure field is most probably wrong, the 

velocities U m and V' will not satisfy the continuity equation. Therefore 


(Sql4) 


the next step is to derive pressure and associated velocity and mass flux 
corrections, so the corrected values will satisfy the continuity and the 
momentum equations. 

If the velocities satisfies the continuity equation, the energy equation can 


P/*GE IS 

POOM^UALiTY 
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be solved in an analogous manner as the momentum equations. 

For turbulent flow, the turbulent kinetic energy * -equation and the 
dissipation rate e -equations are assembled and solved respectively. 

The values of k and e are used to calculate eddy viscosity and hence 

the diffusivides in the finite volume coefficients. 

The residual norms are calculated for all conservation equations and 
normalized by appropriate reference quantities. The convergence criterion 
is checked, and based upon the result, the algorithm returns to (step 1), 


using new values of the dependent variables, or the iteration process is 
terminated. 


(3.3) Code Modification for Oscillating Flow 

To solve all objectives in this work, the following major changes have been made 


to the original C.A.S.T computer program. 

■ The Lam-Bremhorst version of the low Reynolds number k - 1 turbulence 


model was implemented in order to improve fluid flow and heat transfer 
calculations in the near-wall region. Accordingly, changes in the formulation of 
the boundary conditions (in particular the wall boundary), and appropriate 
damping functions used in the k-t equations, are made. 

The oscillatory flow has a very complicated structure, especially in the 
Stirling system, where continuous change from laminar to transition and to 
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^siiiiung or a cycle, 


_ r, iiVlW IUUC UiC 

” “ deSCribed “ ^ ^ -ta «• acceleration of the fluid is 
taJdng place, a turbulent slug forms and advances in the flow direction. The fluid 

“ mbe " 1 Stffl ' * ' *— - — * *h= Portion Of the tube where the 
turbulent dug did no. arrive, and turbulent in the region through which the dug 

passed. The flow is presumed to be turbulent when the leading edge of the 

Slug arrives at a given axial location within the tube. The empirical transition 

mode, used to determine the dug position wili be describe in Chapter VI. 

TO insure iaminar flow in the porrion of the tube where the dug has no, 

yet arrived, the turbuien, viscosity, which is a part of the effective viscosity was 
suppressed (set to zero). 

original C.A.S.T. code (Peric and Scheuerer (1989)), was designed 
to solve steady and unsteady problems. Actually the unsteady case was solved by 
dividing the rime domain into a number of time steps (steady cases), where the 

uutial guess for the following one. Regarding the oscillation, the acceleration and 
^deration of die fluid is accounted for by specifying inflow axial velocity (see 
°3uarion (2.17)). The inflow and outflow boundary planes are switched every 

I8 °“' " the mean vdoci >y from one compile cycle is aero. 
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CHAPTER IV 


Turbulence Modeling of Oscillatory Flow 


Since the early seventies, a number of k-t models have been developed and 

implemented into various engineering applications. In all of these models, the eddy 
viscosity concept is used along with two additional partial differential equations which 
are derived by manipulating the Navier-Stokes equations. These are the transport 
equations for the turbulence energy ( k -equation) , and the isotropic turbulence dissipation 

rate ( e -equation). 

(4.1) Comparison of Various Turbulence Models 
The set of model equations recommended by Launder and Spalding (1972) for 
high Reynolds number flows have been most widely employed. For the wall bounded 
Hows these equations are used in conjunction with the empirical wall function 
formulations. The wall functions translate the wall boundary conditions into the region 


Witte, r - 50 distance from the wall w „ 

0974) eXte " de<1 the*-, 

model to low Reynolds number, and perform^ 

performed computations right to the wall T„ 

other forms of the k-t mnw i aU ‘ Later » 
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the damping functions were #■ 
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of modeIs [--and Lam-Bremhorst] have difficulty when the initial conditions contain 
luge gradients. ..." 

On the basis of these tests, the simplicity of the model, and the successful 
application to the turbulent oscillatory flow problem, Koehler ( 1990 ), the Lam-Bremhorst 
model was chosen for present study. 

(4.2) The Lam-Bremhorst k-z Model 

All of the different forms of k-t models based on the eddy viscosity concept use 

the same form of equations for * and e . The terms of these equations are presented in 
table II. 



Convection 


Diffusion 


.7; dk j»dk 

P U-— + pV — 
dx dr 


+pt/|i + pF~ 

dx dr 


Terms of the k and e equations. 


■ 

= 

dx a k dx 

t or a k dr 

+ pG 

= 

dx a t dx 

r or a t dr 

*c,/,pOj 
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The turbulent viscosity i s modeled • 

’ Usui S Prandtl Kbmogorov’ 


s expression as: 
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model with wall function formulas, suggest that /„ should approximately be 

equal to unity in the fully turbulent region remote from the solid walls. This is 
also consistent with the usual understanding of turbulence, that properties should 
be fairly uniform in regions where viscous effects are small compared to turbulent 
ones. On the other hand in regions very near a wall where viscous effects become 
important properties will change rapidly and /j, will also differ considerably from 


unity. 

Function f x 

Computations with the high Re number form of the model with wall function 
formulas suggest that f x is approximately unity remote from the wall. In the near 


wall region it is founded that /, assumes larger values in order to increase the 


predicted dissipation rate, thereby reducing the predicted turbulence level to 


match available experimental data. Lam and Bremhorst proposed that: 


A 


W (^) 3 


( 4 . 4 ) 


where f x is a function of /j, only, with constants obtained by computer 
optimization. Close to the wall, will be small but finite, and f x will become 
large. Away form the wall the turbulence level and /j, are high. Hence, f x will 
be approximately equal unity. 
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Figure 4.1. The damping functions in the Lam Bremhorst model. 

(4.2.2) The Boundary Conditions in the k-e Model 

The high Reynolds number model is a special case of the low Reynolds number 
model, where all damping functions are set to unity. The real difference between the high 
Reynolds number and the low Reynolds number models is due to the boundary 
conditions. The boundary conditions are presented in table IV. 

Because in some low Reynolds number models the wall values of k and e are 

defined differently, the additional terms appear in the turbulent transport equations. The 

^ Bremhorst model offers the advantage that there are no additional terms. It makes 
this 

model easy to implement on base which is C. A.S.T. code, where the high Reynolds 
number model was originally used. 




(4.3) Evaluation of the Constants in the k-t Model 



There are five empirical constants used in the k-z turbulence modeling. 
c„ was determined from experiments in thin shear layers using the relation: 

.[/Vs. 2 ( 4 -6) 

c “ “ ( k * 

The value in the above equation was measured by Champagne, Hams and 
Corrsin (1970) as ■ 0.09 . 

c, was found stitdying isotropie turbulence for high Re number flows. Batchelor 
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constants used in this work are listed 


in table V. 


Table V. Turbulent modeling 


constants. 
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arrival of this slug to any axial location, the flow becomes turbulent. The following is 
used to find an expression for the point in the cycle at which the leading edge of the slug 
appears at the particular axial position. Using definition: 

= ^sin(0) 

t 

0 

and substituting 0 = <of , Re^ = t/^D/v and Va = coD 2 / 4v gives: 







t 



0 



** V 


£><•> 


(l-cos(or)) 


X *J° - |-^(l-cos(8)) 


where -y.,„ g is an axial slug position. 


CHAPTER V 

Code Vacation avith Pipe How under Steady State Conditions 

In order to verify bod, turbulence modeling and » -purer code, a series of 

computational tests were made to compare the pred.cuons for pipe 

conditions. Concerning fluid mechanics, calculations were made a. four t 
inflow velocities to cover laminar, transitional and turbulent flows, using three numenc 

models. These models arc. 

•tin tVu* fluid where the molecular viscosity 
shifted away from the wall to the pornt m the fluid wne 

effect is small; 

■ die Low Reynolds Number *-« Model, where the turbulence transport 
equations are solved, the wall damping effect is modeled, and the direct influence 
of molecular viscosity is accounted for; 

. die Laminar Model, where there is no turbulence modeling and turbulent 
intensity is assumed to be zero. 
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universal profile line near the wall, and are shifted away from this line away from the 
^•all outside the boundary layer where the core velocity is constant. The High Reynolds 
Humber Model is not applicable for such a slow flow and predicts a velocity as much too 
high. This model was developed for turbulent flows so such a weak agreement can be 

expected. 

Transition flow: Re = 5000 

Figure 5.2 shows dimensionless velocity profiles for transition flow at Re = 5000. The 
Laminar and the Low Reynolds Number Models show laminar profiles where most of 
the velocity distribution foUows the U* = Y* curve. The High Reynolds Number Model 
predicts turbulent flow where the velocity distribution roughly follows the 
U' = 2.44 ]n(Y*) + 5.5 curve. None of the solutions for steady state flow at Re = 5000 
appear to be satisfactory, as seen in figure 5.2. 

Turbulent flow: Re = 15000 

Figure 5.3 shows velocity profiles in dimensionless coordinates at the inflow velocity 
corresponding to Re = 15000 . At this value of the Reynolds number, the flow is 
turbulent. At this condition, the Laminar Model performance is not good because it 
predicts a laminar profile according to the relation U* = Y* not only in the near wall 

region, but also throughout the core flow. The High Reynolds Number Model follows 
the universal profile in the logarithmic region, but and is not correctly predicting the 
velocities in the viscous region where the wall functions are used. The Low Reynolds 
Model is seem to match the universal profile very well, and is the only model which 



correctly works near and away from the wall. 
Highly turb ulent flow! Re = 50000 


Figure 5.4 shows dimensionless velocity profiles for turbulent flow at Re = 50000. The 
comparison between the discussed numerical models is the same as for turbulent flow 
with Re = 15000 . As expected, the Laminar Model shows its inapplicability for 
turbulent flow. The Low Reynolds Model matches the universal profile near the wall, in 
the buffer zone, and in the logarithmic and the overlap regions. The High Reynolds 
Model predicts the velocity profile only in the logarithmic and overlap regions. 

From the above discussion, one can see that the Laminar Model, where only 
molecular viscosity is used and no turbulence is modeled, can be used for flows in the 
laminar and the transition regimes. Only in these regimes, the results are comparable 
with the more advanced U>w Reynolds Model. The High Reynolds Model shows good 
performance in handling turbulent flows, but yields incorrect results in the laminar cases. 
Among tested models only Low Reynolds Model shows an ability obtain accurate 
velocity profiles in steady state laminar, transition, and turbulent flows. 
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Figure 5.3. Dimensionless velocity profile for the steady state 
fully developed pipe flow at Re - 15000 . 
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CHAPTER VI 

Numerical Results and Comparison with Experiment 


Modeling efforts for turbulent fluid flow and heat transfer ean be evaluated by 
comparing the computational predictions with experimental dab,. The code valrdabon 
section (Chapter V) provides the basis of these 'valuations and establishes 
confidence based no, only on the computation scheme used but also u^n die turbulence 
model utilized. In this thesis, turbulence modeling assumptions for steady state conditions 
have been extended to unsteady flow conditions, and particularly oscfllatory (with zero 
mean) flow conditions. In this chapter the numerically predicted solutions are com 
with experimental dam. obtained from an oscillatory flow apparatus at die University of 

Minnesota, which now will be briefly described. 
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10 measure the mean velocity, U, and the rms fluctuation of the axial component of the 
velocity, U ' . The cross-wire probes were used to add the radial velocity components 

mean, V, and rms-fluctuation, V , as well as the Reynolds shear stress, -U'V* . The 

Top Dead Center (TDC) a photodetector was used to detect the position of the piston, 
so that the hot wire measurements could be related to a certain crank angle. Typically 
the measurements taken were an ensemble average of the measurements obtained over 
500 cycles or more. 

The experiment was designed and operated to simulate the SPDE Stirling engine 
heater performance in terms of three dimensionless parameters, the maximum tube 
Reynolds number of the cycle (Re^), dimensionless frequency of oscillation (Valensi 

number Va ), and mean fluid displacement rates (A r ). Air at room temperature was used 
as a working fluid. The operating parameters of the experiment are listed in table VI. 



11840 

Va 

80.2 

K 

1.22 

D [ m ] 

0.038 

I/D 

60 


Table VI. The operating parameters. 

For a more detailed description of the experiment and for all the experimental data refer 
to Seume et al. (1992). 
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(6.2) Fluid Flow Predictions 

Three numerical models were used to predict the fluid flow: 

■ THE LAMINAR MODEL, where the flow is assumed to be laminar everywhere, 
so no turbulence modeling is needed. In the numerical code the value of turbulent 
viscosity was set to zero ( \i T = 0 ). 

■ THE TURBULENT MODEL, where turbulence is modeled utilizing the Lam- 
Bremhorst version of a low Reynolds number k-t turbulence model. In this 
model the k-t is kept active throughout the cycle and at all axial locations. 

■ THE TRANSITION MODEL, where an empirical transition model has been 

utilized to activate the k-t model at different times of the cycle and axial 
locations in the tube. 

In the computation a 60 (axial) by 70 (radial, half of the pipe) grid was used for all 
unsteady cases, with the grid density being high near the wall and sparse away from it 
in the transverse direction, and a uniform grid density in the axial direction. The 
convergence criterion was set as a 0. 1 % of the global residual norms for every dependent 
variable. For each case 120 time steps per cycle were used. In most cases the steady 
oscillating flow conditions, when there is no cycle-to-cycle variation, were reached after 
ihrce to four cycles. Under laminar flow conditions throughout the cycle (not the focus 
°f this study - see Ahn and Ibrahim (1992)) the comparison between experiment and 
computation should be done using the steady oscillating flow results, i.e.in the fourth 
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For a,, c “ P 3 “' RgUre 6 2a ahows similar p,ots to Figure 6.2b at X,D = , 6 

- — — «• — . 1 ' 
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^ment is SCO* except -the legation of — — ■ ~ 

predictions, are good at locations for 8 x 150* , but different from me expert mental^ 

l crank angles 8 - 170* , where the flow reversai near the wail takes p>ace. Figure 
6 . 2C shows sintiiar piots to Figure 6.2b at XI D - 44. The agreement between the 
experimental dam and the computational results are similar to what was desenbed m 

Figure 6.2b. 

Comparing Figures 6.2a,b,c one can see dm. excellent agreement between the 
experimental data and me computarions in dte laminar portion of the cycle takes place. 
Ue occurrence of transition is located accurately by the Transition Model. More over, 
agreement between dte experiment and computarions in the turbulent poruon of dte 

cycle is good, and gets better at higher XID values. 


TWnrinnless Velafltt Profiles (11+ VS X+ 1 

Figures 6.3a, b,c show dte dimensionless velocity V vs dimensionless distance 

from the wall Y' on a semi logarithmic scale for selected crank angles at XI 
and 44 . This form of dam presentation allows looking at the velocity near the wall in 
dm viscous sublayer region. On dte plots, the universal velocity profile from me law of 
me Wall (Iri - r , cr - «*>.r ♦») is p-nmd (dashed line). Comparison 
between me experimental dam and me Transition Model predictions is summarired in 


Table VD. 



Location 

in 

the Tube 

g 

Flow Type 
from 

experiment 

Computational 
Agreement with 
experiment 

X/D = 16 
X/D =30 
X/D =44 

30 

laminar 

laminar 

laminar 

excellent 

excellent 

excellent 

X/D =16 
X/D =30 
X/D =44 

60 

laminar 

laminar 

laminar 

good 

good 

good 

X/D = 16 
X/D =30 
X/D =44 


turbulent 

turbulent 

laminar 

excellent 
very good 
good 

X/D = 16 
X/D =30 
X/D =44 

HI 

laminar 

turbulent 

turbulent 

fair 

very good 

vprv rrivvl 

m 

RRWMl 


gtXXJ 

mmmm ■ 

laminar 

turbulent 

turbulent 

fair 

good 

excellent 


turbulent 

Table VII. Evaluation of Figures 6.3a, b,c. 
e: eXCeIIenl ±5% - vei y *°°d ±10* , good ±25% , fair >±50% . 

Figures 6.3a, b,c provide another way of comparing the experimental data with 
* & th '“ f ‘ 8Ur “ * «gion is expanded and friction velocity 

“ ° f b0th ^d distance. Simiiar conclusions to these 

from Figures 6.2a, b,c are noticed. Again the Transition Model is capable of predicting 

“>e laminar and turbulent parts of the cycle accutately. The only difficulty is for high 
(near 8 - 170* and above for X/D - 30.44 , and 8 = 120* and above for 
X, ° ~ 16) Where the flow decelerates and relaminarization might occur. 
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Figure 6.2a. Normalized veloc: 
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(6.2.2) Turbulent Kinetic Energy 

In the k-t turbulence modeling, the transport equation for the turbulent kinetic 
energy (equation 2.7), is a basis for describing the transport processes in fluid motion. 
The turbulent kinetic energy is defined as: 

k = (l/2)(l/ /2 + V /2 +^ 2 ) 

In the experiment, where the hot-wire anemometer was used, the fluctuating components 
of the velocity ( U' and V*) were measured. In order to obtain k , one should make an 

assumption about W . In this study, W = V , which is appropriate in the turbulent core 
where turbulence is actually isotropic. 

k = (\l2)(U /2 +2V /2 ) (6 * 2) 

Figure 6.4a illustrates the comparison between numerical and experimental predictions 
for turbulent kinetic energy vs distance from the wall for different crank angles at the 
axial location XI D = 16 . In the laminar portion of the cycle, at 6 = 30° , the 
computed values of A:, as well as experimental data, are close to zero. Later, when a 
turbulent slug is advancing into the tube, the values of k increase. The highest values 
of k are obtained at about 6 = 90° , where the agreement with the experimental results 

is the best. At this crank angle, the flow is turbulent. The calculated value of k rises 
quickly from zero at the wall, up to a maximum near the wall, and gradually decreases 
to its lowest values at the center-line of the tube. For larger crank angles, the 
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m 


experimental results show a decrease in the value of * 


as a result of the flow 


deceleration and possible relaminarization. Figures 6.4b and 6.4c show plots similar to 
those in Figure 6.4a but for X,D - 30 and 44 respectively. A, X,D - 30 and 44 
good agreement between predictions and experiment are noted at almost all crank angles 
for the value of t. It should be noted that the profiles shown on Figure 6.4b are for 

6 a 60” , the slug arrives at about 0 = 75* , after which the flow become turbulent. 
Similarly, on Figure 6.4c only profiles for 0 a 90’ are shown since the slug arrives at 


about 0 = 105 a 


(6.2.3) Skin Friction Factor 

Figure 6.5 shows the Transition Model prediction of the Motion coefficient as it 
varies through the first half cycle at three locations of the tube, X/D = 16. 30, 44. Also 
on the figure, the experimental data, the results from the Laminar Model and the 
Turbulent Model are shown. The Motion coefficient, C. , is defined as: 


(l/2)p U? 


* 


(6.3) 


“here V, is the average instantaneous velocity across the cross section of tire channel. 

a, the crank angles 0 = 0-.180” tire average instantaneous velocity values are 
“taost zero, and the Motion coefficient approaches infinity. In the laminar portion of the 
^ole, tire laminar Model predictions are in excellent agreement with tire experimental 
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Transition Model 
- - Experiment^ 







Figure 6.4c. Turbulent kinetic energy (Jfc) profile for oscillatory pipe flow at 

XI D = 44 . 









dau. In the turbulent portion of the cycle the Turbulent Model predicts a sldn 
friction factor which is in good agreement with the experimental data. Using this mode, 
the postdoc of transition is predicted a. the same crank angle in the whole tube, a. about 

9 - . (as seen in Figure 6.5), which is no, in agreement with the experiment 

Utilizing the Transition Model, the location of transition is predicted accurately in the 
tube. Moreover, the agreement between the model's prediction for C/ and the data is 
excellent in the laminar par. of cycle and gcrnd in the turbulent part of the cycle. Figure 

6 6 Sh ° WS Pl0 “ ^ 10 6 ' 5 ’ «*» <* friction velocity U' is presented 

instead of the skin friction factor. By definition 


u m = 




li 

P 


(6.4) 


the friction velocity is proportional to the 


square root of the wall shear stress value, 


" ^ ^ ^ ^ as the normalisation factor as in the 

definidon of the skin fricdon factor (e,. 6.3). Therefore values of W are Cose to zero 

« crank angles near 6 - 0-.U0- . Moreover, the agreement betwemt the Transidon 
Model and the experiment is excellent no, only in the laminar portion of cycle bu, also 


in the turbulent portion. 


From die results presented on figures 6.5 and 6.6 and from the above discussion 
cn see tha, the numerically predicted values rue in excellent agreement with the 

experiment. Once agam, the Transidon Model shows the ability to accurately determine 
the Iocado " °f transition throughout the cycle. 


Skin Friction Factor, Cf 



Figure 6.5. Skin friction factor ( C,) predictions front the Turbulent Model, 
the Transition Model, the Laminar Model, 
and the experiment for oscillatory flow 
(Re ^ = 11840 , Va = 80.2 , LID = 60) at 
X}D = 16 , 30 , 44 locations. 






Friction Velocity, u* 



Figure 6. 6. Friction velocity (V \ 
^ Transition 


the Ti^stiion Model fr ° m ^ 

and the exDerimen!’^ Model, 

<*u - ns^r/r^ flow 

W.M 30 J; £/ °-«»a t 

* J0 » 44 locations. 


the Turbulent M 

irliil 




(6.3) Heat transfer 

(6.3.1) Description of the Experiment 

The experimental setup at the University of Minnesota has been utilized to obtain 
heat transfer data. In addition to the fluid flow experimental setup described in Section 
(6. 1), two heat exchangers were installed at the ends of the tube with a heating element 
wrapped along the tube in between. The temperature measurements were made at three 

different axial locations X/D = 1 , 11 and 31 , along the tube length ( L = 62 D). Via 
accurate temperature measurements made near the wall at all crank angles, and by 
extrapolating those temperature profiles, the wall temperature and the wall heat flux were 
obtained, (see Simon (1993) for further details). The wall temperature values are 
summarized here in Table 6.2. Also, inlet working fluid temperature was measured and 
found to be constant throughout cycle ( T wflgw = 21°C). 


Axial 

Position 

X/D = 1 

X/D = 11 

X/D = 31 

Wall 

Temperature 

28° C 

36° C 

40° C 


Table 6.2. The measurements of the wall temperature. 


(6.3.2) Numerical Predictions 

From several computational experiments, it was found that the numerical 
results are very sensitive to the temperature at the wall boundary, and to the initial 





temperature inside the tube. From experimental data obtained from (Simon (1993)) 

and listed in Table 6.2) two different estimates of wall temperature distribution, 
presented on Figure 6.7, were chosen. 



WALL BOUNDARY CONDITIONS (BCI) 



WALL BOUNDARY CONDITIONS ( BC H ) 


Figure 6.7. Wall temperature distribution. 

•n« boundary condition (BC I), is defined by connecting the experimental temperature 
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points with straight lines, assuming the maximum value of wall temperature 
( ^wau = 40 ) occurs at the midpoint of the tube length. On the other hand, the boundary 

condition (BC II), is assumed to experience a temperature plateau, with T^ u = 40° in 

the middle one/third of the tube. In regard to the initial temperature profile inside the 
tube, the situation is more complex than in the fluid flow case. The fluid flow 
calculations were conducted assuming that complete mixing was taking place at the end 
of the cycle and a uniform inflow velocity distribution was adequate. Accordingly, the 
laminar flow calculations were not affected by the previous half cycle turbulent 
calculations. However, the temperature field calculations at the end of first cycle depend 
upon the amount of turbulence left in the tube at that time. In the next half cycle, that 
temperature field taken as the initial condition, affects the laminar (and later on the 
turbulent) thermal field. Therefore, several cycles are necessary to obtain a complete 
thermal field solution with the understanding that the laminar flow thermal field is 
affected by the previous half-cycle turbulent thermal field. 

(6.3.3) Temperature Profiles 

Figure 6.8 shows the midplane temperature profiles from computations and the 
experimental data at different crank angles, using BC I in this computation. (Only a slight 
difference was noticed in the temperature profile using either BC I or BC n.) The 
predicted values of temperature match the experiment close to the wall, but in the core 
flow these values are lower than the experimental data. Near the wall the predictions are 



in close agreement with the experiment since the computed wall temperatures match the 
values of wall temperature obtained experimentally. However, in the core of the tube the 
disagreement is attributed to the turbulence model used. Notice that the turbulent 
temperature predictions (e.g. 0 = 90° ) are not in agreement with the data. 

Consequently, the end of the cycle thermal field (also the initial thermal field of the 
following cycle) does not match the experiment. This results in disagreement between the 
computational results and the experimental data in the laminar portion the of cycle. 

In addition the temperature profiles at two other axial locations are presented in 
Figures 6.9a and 6.9b. Because the relative amplitude of the fluid displacement is 
(A r = 1.22 > 1.0 ), which means that the fluid is displaced more then the length of the 

tube, the slug of cold fluid penetrates most of the tube. At XI D = 16 , the temperature 
at the center line drops down to the inflow temperature = 21° C) after 0 = 30 . 

At XI D = 44 , in the laminar portion of cycle, the centerline temperature rises, later 
drops down when a cold turbulent slug arrives, and rises again when the turbulent flow 
decelerates. The numerical results in the laminar portion of cycle can be further 
improved by making the temperature field at the beginning of the cycle more accurately 
reflect the actual running conditions in the experiment. 

(6.3.4) Wall Heat Flux 

Figure 6. 10 shows the calculated and experimental data of the wall heat flux vs 
crank angle at the midplane. The calculations are made with the Transition Model using 
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two different boundary conditions (BC I and BC II) generated from wall temperature 
measurements. The agreement with the experimental data in both cases is excellent in the 
turbulent flow regime and good in the laminar flow one. Using BC I, higher values of 
the wall heat flux were obtained. This is expected, because BC I implies that fluid arrives 
at the midplane with a lower temperature than when using BC II. At the beginning and 
end of each half cycle the difference in the heat flux between BC I and BC II is small. 
Figure 6. 1 la shows the calculated values of the wall heat flux vs crank angle, at the axial 
location XI D = 16 . At this location the flow becomes turbulent at a lower crank angle 
(i.e. earlier), so the maximum values of wall heat flux for both boundary conditions are 
shifted to lower crank angles. Figure 6. 1 lb shows the wall heat flux data vs crank angle 
at XI D = 44 . At this location the heat flux values are lower, because the arriving fluid 
is warmer than at the midplane. This low heat flux occurs despite the fact that the wall 
temperature is lower than T nvll = 40° which occurs, (see Figure 6.7) at the midplane. 

Using experimental data documented in Seume et al (1992), and Simon and Qiu 
(1993 a,b), a one dimensional model (1-D Model) for estimating laminar and turbulent 
wall heat flux in oscillatory pipe flow has been formulated by Ibrahim et al (1993). In 
this model for the laminar portion of the cycle, the wall heat flux is calculated using the 
Smith and Spalding method as: 



n< *-(T w -T a .) 
1.43 Pr Va 0J - — 

D 


sm 3 (0) 


\ 2 + cos 3 (0) - 3cos(0) 


(6.5) 


For the turbulent heat flux, the largest value form the following two equations has been 


used: 


It = 0.Q21 Pr 06 A (7 > ^ . , Q} 


( 6 . 6 ) 


q ” = 003 P c ^(^-^ r ) (6.7) 

where (J^) and (7^), are the laminar and the turbulent sink temperatures, (see 

Ibrahtm et. al. (1993)). Figure 6.12 shows the 1-D Model predictions for the wall heat 
flu, a. the midplane of the tube. In the figure, the experimental data and results from 
““ Transition Mode, are presented. The 1-D Mode, predictions are seen to be good in 
the turbulent pan, while agreement is fair in laminar flow. As expected the Transition 
Mode, performance is better than the 1-D Mode,. This 1-D modeOng analysis was done 
.n an attempt to avoid using extensive two-dimensional computations, and with the intent 

masks many physical details that are better described by a 2-D Mode,. However it is 
surprising that the ,-D Mode, solution shows considerable success. 
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CHAPTER VII 


Closure 


( 7 . 1 ) Summary and Conclusions 


1) Numerical scheme 

The code used in .his s.udy is based on ft. C.A.S.T. (Computer Aided 
Simulation of Turbulent Flows) code developed b, Peric and Scheuerer (1989), 
which was modified to solve for unsteady oscillatory flow with zero mean 
velocity. Regarding the oscillation, tire acceleration and deceleration of the fluid 
is accounted for by specifying tit. inflow axial velocity, and the flow in the inlet 

and the outlet boundary planes is switched a. the half cycle (180‘) and the end 

of each cycle (360°). 

(2) Tnrhnlence modeling 

The low Reynolds Humber (lam-Bremhors.) k-t turbulent model was 


introduced, in order to improve fluid flow and heat transfer calculations in the 

near wall region. The High Reynolds number *- e model, which was originally 

used m the C.A.S.T. code was found not be accurate enough for solving the 

laminar and turbulent portions of the oscillatory flow. Moreover, an empirical 

transition model was utilised to activate the low Reynolds turbulence model at the 

appropriate time within the cycle for the axial location within fee tube. In fee 

computations a 60 (axial) by 70 (radial, half of fee pipe) grid was used for all 

unsteady cases with fee grid density being higher near the wall and lower toward 

the core, and a uniform grid density in fee axial direction. The convergence 

criterion was se, as a 0.1% of the global residual norms for every dependent 

variable. For each case 120 time steps per cycle were used. A typical run 

involved approximately 2000 seconds of CPU time per cycle on a Cray XMP 
supercomputer. 

(3) Code validation 


The developed numerical code performance as well as the turbulence model was 
validated for steady-state fluid flow and heat tiansfer cases. The study was 
conducted for steady state pipe flow wife uniform inflow velocity distribution for 

four values of Reynolds number (Re - 500 , 5000, 15000, and 50000 ). 

To insure feat felly developed conditions were achieved, a tube length of 120 
diameters was chosen. The low Reynolds Number turbulence model was found 
•o he in excellent agreement wife fee universal velocity pmfifc from the law of 

the f0r *“ fl0W types studied - laminar Model (wife no turbulence 
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modeling), and the High Reynolds Nuraber ^ ^ ^ ^ 

perfonn well in laminar and turbulent flow types respectively. In the steady state 

flow ** - — - - — >-„ „„ m , r va j were 

compared with the analytical results. The smdy war conducted for steady «, 
laminar flow (* - 500), and transition flow (A = 150 00) fa die tube with 
constant wall *mpeiature. To insu re a fiiUy developed temperature profile, the 
■ube with a length 300 diameters was chosen. The value of the Nussel, number 
predicted by the laminar Mode, for die laminar flow case was greater by6% 
•ban the analytical solution. The low Reynolds *-e mode, showed die same 

accuracy in calculadng Mr for die laminar flow case. For the turbu,ent flow case 
difference between die Low Reynolds Number Turbuknce ^ ^ ^ 

analytical solution was 20% . Such a Urge difference can be attributed to 
inaccuracy of the turbulence modeling. 

( 4 ) Unsteady fln.vi 

TTc oscillating flow dimension, ess parameters, ^ 11840> ^ ^ 

, LID - 60, and d, - 1.22 , were chosen to match the Space Power 

Demonstration Engine (SPDE) operating conditions. The computational results 

for: normalized velocity profiles (Ul(J g vs r/R), (2) dimensions 

velocity profile ( U* vs Y* \ n\ ^ . 

), (3) turbulent kinetic energy (k vs r /J?), ( 4 ) 

skin friction factor (C a ^ „ . 

/ )>and (5) friction velocity ( U' vs 6), were 
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compared with the experimental data at three axial locations within the tube 
(Xf D = 16, 30 , and 44). The agreement between the Low Reynolds Number k-t 
model, with an empirical transition model, and the experiment is excellent in the 
laminar portion of the cycle and good in the turbulent portion. Moreover, the 
location of transition has been predicted accurately. 

(5) Unsteady fluid flow with heat transfer 

The heat transfer predictions were made for the same fluid flow conditions as 
those described above, with constant fluid temperature at the tube inlet and 
nonuniform temperature distribution at the tube wall. The predicted temperature 
profile and wall heat flux at the midplane were compared with the experimental 
data. Near the wall, the predicted temperature profile was in close agreement 
with the experimental values, but in the core of the tube some disagreement was 
noticed due to difficulty in defining the initial and boundary conditions of the 
fluid temperature and to the turbulence model used. The wall heat flux 
calculations were found to be in excellent agreement with the experimental data 
in the turbulent flow regime and in good agreement in the laminar one. 

(6) Comparison with 1- D model 

Since the wall heat flux predictions match the experimental data, the Low 
Reynolds Number k-t model, with the empirical transition model, can be used 

for testing the much simpler and less accurate one-dimensional models (used for 
1-D Stirling Engine design codes), by generates wall heat flux values at different 
operating parameters than those of the experimental conditions used herein. 


81 


\ ® ecommen dations for Future Research 

The principal objective of this researrh 

of predicting heat transf * CVel ° P a nume ncal model capable 

exchanger some suggested area? f 5 engine heat 

jested areas for ftture studies are outlined below; 


□ 


□ 


JT" “■ ta — * — — , « „ 

K ° CC “ rS " ““ " ,d ° f rach ** cycle can be modeled 

the SUrlln e engine such as the- heater r 

■ neater ' regenerator and cooler. 


One can always expect to find other areas nf - , 
author anticipates that the data provided fro thi "" ^ "* Pr ° b ' em ' ^ 
investigations on this sn,ect. ~ ^ ^ 
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